Emergence and structure of decentralised trade networks around dark web marketplaces

Dark web marketplaces (DWMs) are online platforms that facilitate illicit trade among millions of users generating billions of dollars in annual revenue. Recently, two interview-based studies have suggested that DWMs may also promote the emergence of direct user-to-user (U2U) trading relationships. Here, we carefully investigate and quantify the scale of U2U trading around DWMs by analysing 31 million Bitcoin transactions among users of 40 DWMs between June 2011 and Jan 2021. We find that half of the DWM users trade through U2U pairs generating a total trading volume greater than DWMs themselves. We then show that hundreds of thousands of DWM users form stable trading pairs that are persistent over time. Users in such stable pairs turn out to be the ones with the largest trading volume on DWMs. Then, we show that new U2U pairs often form while both users are active on the same DWM, suggesting the marketplace may serve as a catalyst for new direct trading relationships. Finally, we reveal that stable U2U pairs tend to survive DWM closures and that they were not affected by COVID-19, indicating that their trading activity is resilient to external shocks. Our work unveils sophisticated patterns of trade emerging in the dark web and highlights the importance of investigating user behaviour beyond the immediate buyer-seller network on a single marketplace.


S2 Dark web marketplaces and identification of real identities performing Bitcoin transac-
tions.

S2.1 Dark web marketplaces
DWMs are in many ways similar to other online marketplaces. They have strict policies that every user must follow. For instance, in some DWMs are banned categories of products, like human trafficking, contract killing, weapons, or COVID-19 fake vaccines 8,9 . Registration is required for all sellers, and sometimes also for buyers. Certified sellers can advertise their products. They have a reputation, which is based on buyers' reviews 10,11 . They are also responsible for delivering the products, sometimes with a tracking number attached, and may offer refunds or reshipment. Buyers are free to look at the listings and sometimes can ask questions directly to the relative seller 12,13 . Payments are often protected by escrow services. These are third-party services, which guarantee that buyers can safely have their money refunded. Users' on DWMs constitute an active community. Numerous are websites and forums where users can share their experience and get advice on the most trustworthy DWMs and sellers, such as Dread 14 , Raptor.life 15 , DarkNetLive 16 , and DarkFail 17 .
DWMs have some unique features as well. They sell several kinds of illicit products, like drugs, fake IDs, and medicines [18][19][20] . They are not accessible by standard web search-engines, but operate online in an encrypted part of the Internet 21 . Potential buyers can easily access to DWMs using specialized browsers, like The Onion Router (Tor) 22 , and anonymously trade illicit goods using cryptocurrencies, like Bitcoin 23 . Bitcoin is currently the most popular cryptocurrency on DWMs [24][25][26] and its adoption is growing in the regular economy as well. Its infrastructure seems to ensure complete anonymity to its users. If a proper technique is adopted, however, there are chances to link the Bitcoin blockchain (that is, the entire Bitcoin transaction history) with the user's real identity 3 . When the Bitcoin blockchain is successfully linked to a real identity, the records of past, present, and future Bitcoin transactions is traceable, easily accessible, and can be used by companies, law enforcement agencies, and researchers.

S2.2 Identification of real identities performing Bitcoin transactions
The raw, anonymized Bitcoin blockchain can be publicly accessed through Bitcoin core 27 or third-party APIs such as Blockchain.com 28 . It contains information about origin and destination addresses, as well as time and amount of the transactions. In order to contrast traceability of the real identity, an user is likely to use multiple addresses. A new address is often generated in each transaction. Grouping the addresses in clusters reduces the complexity of the Bitcoin blockchain and challenge users' anonymity 29 . Given that millions of Bitcoin addresses are currently active and many others are continuously being generated, a clustering approach primarily based on manual annotation is not feasible. Various heuristics, instead, have been proposed [29][30][31][32] . They were successful in grouping Bitcoin addresses and associate them to cluster of real entities. For instance, in 29 , the authors were able to find a connection between a set of large transactions and a single one, which was dated in November 2010. In 30 , the authors applied to a daily university setting the privacy protocol recommended in Bitcoin transactions, finding that almost 40% of the real identities would be recovered. Another work showed the presence of "super clusters" of entities, which marked macro-variations in the evolution of the Bitcoin economy 31 . The primary reasons behind the effectiveness of heuristic clustering are: "address reuse, avoidable merging, super-clusters with high centrality, and the incremental growth of address clusters" 32 . Figure S1. Identification of real entities in the Blockchain. End goal of Bitcoin transactions clustering techniques: mapping a series of Bitcoin addresses to real entities. In this example, an address sends Bitcoins to another address. Thanks to the identification process, the two addresses are associated with two real entities. The Bitcoin transaction between the two entities becomes traceable and transparent.
The end goal of clustering Bitcoin addresses is to map them to single, real entities, as shown in Figure S1. To achieve this goal, however, heuristic clustering techniques should be improved. Manual annotation has shown a valuable potential 33 . It consists on gathering publicly available Bitcoin addresses, like the Wikimedia Foundation one 34 , and engage through direct interaction with unknown Bitcoin addresses. If some real entities are known, it is easier to associate the remaining Bitcoin addresses to other real identities. In the last few years, companies specialising in Bitcoin analytics have started to leverage previous methodologies [29][30][31][32][33] to unveil real entities. The leading company in analysing Bitcoin transactions on DWMs is Chainalysis Inc. 35 , which has also aided several federal investigations. For instance, it supported the United States Internal Revenue Service (IRS) in tracking Bitcoin transactions 36 and the FBI in the Twitter hack 37 . Chainalsysis clusters Bitcoin transactions in groups by combining previous methodologies [29][30][31][32] and real entities are unveiled with an approach similar to 33 . In the dataset, real entities represent DWMs, users of DWMs, or other real entities interacting with these users. Chainalysis aims at minimizing the false positives, who may lead to wrongly associate a real entity with illicit activities. If a Bitcoin address cannot be uniquely ascribed to a real entity, it is included in our dataset as an independent and unnamed entity. Only a fraction of the entities in our dataset thus represent named and real entities, which identity is known. Given that there are millions of entities in our dataset, it is impossible to identify all the corresponding real identities. After the identification process is completed, to each real entity is associated a string of numbers and the dataset re-anonymized. Transactions to and from Bitcoin trading exchanges are also removed, because our primary interest entails the study of direct interactions between real entities. As a result, the dataset analysed in this article is created by using state-of-the-art clustering techniques to identify addresses owned by the same user and label real entities.  Active in 2020-2021 Closed before 2020

S4 Detection of stable pairs in temporal and directed networks
Here, we summarize the metholodogy of detecting the backbone of stable pairs in temporal and undirected networks as introduced in 38 , and show how it can be easily adapted to tackle the analysis of directed temporal networks. The methodology follow three sequential steps: (i) determine the interval partition, (ii) estimate models' parameters, over successive intervals, and (iii) run a statistical filter, which removes all pairs explained by the null hypothesis and retain stable pairs. The analysed temporal network, either directed or undirected, of N nodes evolves in an observation window composed of T ≫ 1 time steps, labeled as t = 1, ..., T . At each time step t, entities interact among themselves and form a time-varying network of interactions, described by a binary adjacency matrix that varies in time A(t).

S4.1 Temporal and undirected networks
Interval partition. The overall observation window is divided in successive and disjoint intervals using an auxiliary method, namely, the Bayesian Block method 39 . It takes as input the total number of temporal pairs created in the entire network at time t where the superscript "ts" indicates that these variables are estimated from the time series and A ts i j (t) is the i jth entry of the estimated adjacency matrix at time t. The Bayesian Block method returns the interval partition, which divides the overall time window T into I disjoint intervals indexed by ∆ = 1, . . . , I, that contain a uniform total number of connections. From the knowledge of the interval partition, the length, τ(∆), of the generic ∆th interval is obtained with the following closure relation: ∑ I ∆=1 τ(∆) = T . Parameter estimation. According to the null hypothesis, pair of entities i and j are expected to interact proportional to the their individual activities at time t. That is, the probability that entities i and j interact at time t is a binomial random variable defined as where a i (t) and a j (t) are piece-wise constant activities, which represent the propensity of creating interactions at time t. The estimation of piece-wise constant activities is carried out analysing each of the I intervals separately. The activity of entity i at time t ∈ t in(∆) ,t in(∆) + τ(∆) − 1 is computed through the following frequency count: where s ts i (∆) and W ts (∆) ≫ 1 are the total number of pairs generated by entity i in the ∆th and the total number of temporal pairs generated in the network in the ∆th interval, respectively. These variables are computed from the adjacency matrix A ts (t), as, s ts appropriate approximation for long time series. The probability that the observed weight, w ts i j , could be explained by the relative expected weight, E [w i j ] in Eq. (4), is computed according to the cumulative function of the Poisson distribution where P (x; E [w i j ]) indicates the Poisson distribution with random variable x and expected value E [w i j ]. Equation (5) represents the p-value α i j : when the p-value is below a pre-defined threshold, the pair i j is significant and included in the backbone network. The same statistical test is repeated for all pairs of entities i j observed at least once in the overall temporal evolution.

S4.2 Temporal and directed networks
With little modifications, the above methodology can be used to filter temporal and directed networks.
Interval partition. The interval partition is obtained by using the Bayesian Block method as above. The total number of temporal pairs created in the entire network at time t is where not pairs are directed, while in Eq. (1) undirected, thereby explaining the different ranges in the summations.
Parameter estimation. In directed networks, the probability that entity i contacts at random entity j at time t is defined as where a i (t) is the activity of entity i at time where s ts out,i (∆), s ts in,i (∆), and W ts (∆) ≫ 1, are the total incoming strength of entity i in the ∆th interval, outgoing strength of entity i in the ∆th interval, and the total number of directed, temporal pairs generated in the network in the ∆th interval, respectively. These variables are computed from the adjacency matrix A ts (t), that is, s ts out, A ts i j (t), and W ts (∆) = ∑ N i=1 s ts out,i (∆). Once the activity and attractiveness are estimated according with Eq. (8), the probability in Eq. (7) can be evaluated. .

8/16
The probability that the observed weight, w ts i→ j , is explained by the expected weight, E [w i→ j ] in Eq. (9), is computed according to the cumulative function of the Poisson distribution Equation (10) represents the p-value α i→ j , which is used to assess whether the directed pair i → j is significant. The same statistical test has to be repeated for directed pairs observed at least once in the overall temporal evolution. For undirected networks, Eq. (10) is equivalent to Eq. (5).